{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from torch.optim import SGD, Adam\n",
    "from pytorch3d.transforms import Rotate, Translate,euler_angles_to_matrix\n",
    "import torch\n",
    "from torch import nn\n",
    "\n",
    "\n",
    "def euler_to_rotation_matrix_zyx(angle_z, angle_y, angle_x):\n",
    "    # Extract the individual angles\n",
    "    #angle_z, angle_y, angle_x = euler_angles\n",
    "\n",
    "    # Compute trigonometric functions for the angles\n",
    "    cos_x, sin_x = torch.cos(angle_x), torch.sin(angle_x)\n",
    "    cos_y, sin_y = torch.cos(angle_y), torch.sin(angle_y)\n",
    "    cos_z, sin_z = torch.cos(angle_z), torch.sin(angle_z)\n",
    "\n",
    "    # Compute the elements of the rotation matrix\n",
    "    R00 = cos_z * cos_y\n",
    "    R01 = cos_z * sin_y * sin_x - sin_z * cos_x\n",
    "    R02 = cos_z * sin_y * cos_x + sin_z * sin_x\n",
    "    R10 = sin_z * cos_y\n",
    "    R11 = sin_z * sin_y * sin_x + cos_z * cos_x\n",
    "    R12 = sin_z * sin_y * cos_x - cos_z * sin_x\n",
    "    R20 = -sin_y\n",
    "    R21 = cos_y * sin_x\n",
    "    R22 = cos_y * cos_x\n",
    "\n",
    "    # Construct the rotation matrix R from the elements\n",
    "    rotation_matrix = torch.stack([torch.stack([R00, R01, R02]),\n",
    "                                   torch.stack([R10, R11, R12]),\n",
    "                                   torch.stack([R20, R21, R22])])\n",
    "\n",
    "    return rotation_matrix\n",
    "\n",
    "class PolicyNetAtt(nn.Module):\n",
    "\n",
    "    def __init__(self,\n",
    "                 input_dim: int,\n",
    "                 policy_dim: int = 2):\n",
    "\n",
    "        super(PolicyNetAtt, self).__init__()\n",
    "\n",
    "        self.num_landmark = int((input_dim - 4) / 4)\n",
    "        print(\"input dim:\")\n",
    "        print(input_dim)\n",
    "\n",
    "        self.agent_pos_fc1_pi = nn.Linear(4, 32)\n",
    "        self.agent_pos_fc2_pi = nn.Linear(32, 32)\n",
    "        self.landmark_fc1_pi = nn.Linear(3, 64)\n",
    "        self.landmark_fc2_pi = nn.Linear(64, 32)\n",
    "        self.info_fc1_pi = nn.Linear(64, 64)\n",
    "        self.action_fc1_pi = nn.Linear(64, 64)\n",
    "        self.action_fc2_pi = nn.Linear(64, policy_dim)\n",
    "\n",
    "        self.relu = nn.ReLU()\n",
    "        self.tanh = nn.Tanh()\n",
    "        self.softmax = nn.Softmax(dim=2)\n",
    "        \n",
    "    def forward(self, observation: torch.Tensor) -> torch.Tensor:\n",
    "        if len(observation.size()) == 1:\n",
    "            observation = observation[None, :]\n",
    "\n",
    "        # compute the policy\n",
    "        # embeddings of agent's position  \n",
    "        agent_pos_embedding = self.relu(self.agent_pos_fc1_pi(observation[:, :4]))\n",
    "        agent_pos_embedding = self.relu(self.agent_pos_fc2_pi(agent_pos_embedding))\n",
    "\n",
    "        # embeddings of landmarkss\n",
    "        estimated_landmark_pos = observation[:, 4: 4 + 3 * self.num_landmark]\n",
    "        print(\"estimated_landmark_pos.shape:\")\n",
    "        print(estimated_landmark_pos.shape)\n",
    "        #landmark_info = torch.cat((estimated_landmark_pos.reshape(observation.size()[0], self.num_landmark, 2),\n",
    "        #                        info_vector.reshape(observation.size()[0], self.num_landmark, 2)), 2)\n",
    "        \n",
    "        landmark_reshape=estimated_landmark_pos.reshape(observation.size()[0], self.num_landmark, 3)\n",
    "        #landmark_embedding = self.relu(self.landmark_fc1_pi(landmark_info))\n",
    "\n",
    "        landmark_embedding = self.relu(self.landmark_fc1_pi(landmark_reshape))\n",
    "        landmark_embedding = self.relu(self.landmark_fc2_pi(landmark_embedding))\n",
    "\n",
    "        # attention\n",
    "        landmark_embedding_tr = torch.transpose(landmark_embedding, 1, 2)\n",
    "\n",
    "        # mask\n",
    "\n",
    "        mask = observation[:, - self.num_landmark:].unsqueeze(1)\n",
    "        print(\"mask.shape:\")\n",
    "        print(mask.shape)\n",
    "        attention = torch.matmul(agent_pos_embedding.unsqueeze(1), landmark_embedding_tr) / 4\n",
    "        attention = attention.masked_fill(mask == 0, -1e10)\n",
    "\n",
    "        att = self.softmax(attention)\n",
    "        landmark_embedding_att = self.relu((torch.matmul(att, torch.transpose(landmark_embedding_tr, 1, 2)).squeeze(1)))\n",
    "\n",
    "        info_embedding = self.relu(self.info_fc1_pi(torch.cat((agent_pos_embedding, landmark_embedding_att), 1)))\n",
    "        action = self.tanh(self.action_fc1_pi(info_embedding))\n",
    "        action = self.tanh(self.action_fc2_pi(action))\n",
    "\n",
    "        if action.size()[0] == 1:\n",
    "            action = action.flatten()\n",
    "\n",
    "        #scaled_action = torch.hstack(((1 + action[0]) * 2.0, action[1] * torch.pi/3))\n",
    "        scaled_action=action\n",
    "\n",
    "        return scaled_action\n",
    "\n",
    "\n",
    "class Agent:\n",
    "    \n",
    "    def __init__(self,s=torch.zeros([4]),p_b=torch.zeros([3]),p_T=torch.ones([3]),max_num_landmarks=128,fov=torch.pi/2,lr=1e-6):\n",
    "        self._s=s\n",
    "        self._p_b=p_b\n",
    "        self._p_T=p_T\n",
    "        self._max_num_landmarks=max_num_landmarks\n",
    "        self._num_landmarks=torch.randint(low=1, high=self._max_num_landmarks, size=(1,)).item()\n",
    "        self._p_f_W=torch.ones([3,self._num_landmarks]) #feature position in world frame\n",
    "        \n",
    "        self._V_flatten=torch.tensor([100.,1.,1.])\n",
    "        self._V=torch.diag(self._V_flatten)\n",
    "        self._V_inv=torch.diag(1/self._V_flatten)\n",
    "        \n",
    "        #self._R_W_C=zy_rotation_matrix(angle_z=self._s[0],angle_y=self._s[1]) # camera pose in World frame\n",
    "        self._input_dim=max_num_landmarks*4+3+3+4\n",
    "        self._fov=fov\n",
    "        \n",
    "        self._policy=PolicyNetAtt(input_dim=self._input_dim)\n",
    "        self._policy_optimizer=Adam(self._policy.parameters(), lr=lr)\n",
    "        self.loss=0.\n",
    "    \n",
    "    def update_rotation_matrix(self):\n",
    "        self._R_W_C=euler_to_rotation_matrix_zyx(self._s[0],self._s[1],torch.tensor(0.0))\n",
    "    \n",
    "    \n",
    "    def calculate_visibility(self,a,bearing):\n",
    "        a_vis=torch.tensor(a)\n",
    "        v=1/(1+torch.exp(-a_vis*torch.cos(bearing*torch.pi/self._fov)))#*(1+torch.exp(-torch.tensor(a_vis)))\n",
    "        v_max=1/(1+torch.exp(-a_vis))\n",
    "        v_min=1/(1+torch.exp(a_vis))\n",
    "        \n",
    "        v_normalized=(v-v_min)/(v_max-v_min)\n",
    "        return v_normalized\n",
    "        \n",
    "    def calculate_target_visiblity(self):\n",
    "        p_T_C=self._R_W_C.T @ self._p_T\n",
    "        bear_T=torch.acos(p_T_C[0]/torch.norm(p_T_C))\n",
    "        vis_T=self.calculate_visibility(1.,bear_T)\n",
    "        \n",
    "        return vis_T\n",
    "    \n",
    "    def calculate_FIM(self):\n",
    "        p_f_C=self._R_W_C.T @ self._p_f_W\n",
    "        \n",
    "        M_total=torch.zeros([3,3])\n",
    "        for i in range(self._num_landmarks):\n",
    "            v_f_C=p_f_C[:,i]\n",
    "            dis_f=torch.norm(v_f_C)\n",
    "            x_f=v_f_C[0]\n",
    "            bear_f=torch.acos(x_f/dis_f)\n",
    "            vis_f=self.calculate_visibility(44.,bear_f)\n",
    "            #print(v)\n",
    "            unitX=torch.tensor([1.,0.,0.])\n",
    "            R=rodrigues_rotation(unitX,v_f_C)\n",
    "            \n",
    "            M=vis_f*(R.T@self._V_inv@R)/dis_f/dis_f\n",
    "            M_total+=M\n",
    "            \n",
    "        return M_total\n",
    "    \n",
    "    def plan(self,s,p_b,p_T,p_f):\n",
    "        padding=torch.zeros(3*(self._max_num_landmarks-p_f.shape[1]))\n",
    "\n",
    "        \n",
    "        mask=torch.tensor([True]*(p_f.shape[1]+2)+[False]*(self._max_num_landmarks-p_f.shape[1]))\n",
    "        net_input=torch.hstack([s,p_b,p_T,p_f.flatten(),padding,mask])\n",
    "        action=self._policy.forward(net_input)\n",
    "        \n",
    "        print(mask)\n",
    "        return action       "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Tensor created on CUDA: tensor([1, 2, 3, 4, 5, 6, 7], device='cuda:0')\n"
     ]
    }
   ],
   "source": [
    "# Check if CUDA is available\n",
    "if torch.cuda.is_available():\n",
    "    device = torch.device(\"cuda\")  # CUDA device object\n",
    "    tensor = torch.tensor([1, 2, 3,4,5,6,7]).to(device)  # Creating a PyTorch tensor on CUDA\n",
    "    print(\"Tensor created on CUDA:\", tensor)\n",
    "else:\n",
    "    print(\"CUDA is not available. Tensor will be created on CPU.\")\n",
    "    tensor = torch.tensor([1, 2, 3])  # Creating a PyTorch tensor on CPU\n",
    "    print(\"Tensor created on CPU:\", tensor)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "__init__() got an unexpected keyword argument 'device'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[8], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m ag\u001b[39m=\u001b[39mAgent(device\u001b[39m=\u001b[39;49mdevice)\n\u001b[1;32m      2\u001b[0m ag\u001b[39m.\u001b[39m_policy\u001b[39m.\u001b[39mtrain()\n\u001b[1;32m      3\u001b[0m action\u001b[39m=\u001b[39mag\u001b[39m.\u001b[39mplan(ag\u001b[39m.\u001b[39m_s,ag\u001b[39m.\u001b[39m_p_b,ag\u001b[39m.\u001b[39m_p_T,ag\u001b[39m.\u001b[39m_p_f_W)\n",
      "\u001b[0;31mTypeError\u001b[0m: __init__() got an unexpected keyword argument 'device'"
     ]
    }
   ],
   "source": [
    "ag=Agent(device=device)\n",
    "ag._policy.train()\n",
    "action=ag.plan(ag._s,ag._p_b,ag._p_T,ag._p_f_W)\n",
    "torch.norm(action).backward()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Agent:\n",
    "    \n",
    "    def __init__(self,device,s=torch.zeros([4]),p_b=torch.zeros([3]),p_T=torch.ones([3]),max_num_landmarks=67,fov=torch.pi/2,lr=1e-6):\n",
    "        self._device=device\n",
    "        self._s=s.to(device)\n",
    "        self._p_b=p_b.to(device)*2\n",
    "        self._p_T=p_T.to(device)*3\n",
    "        self._max_num_landmarks=max_num_landmarks\n",
    "        self._num_landmarks=torch.randint(low=1, high=self._max_num_landmarks, size=(1,)).item()\n",
    "        self._p_f_W=torch.ones([3,self._num_landmarks]).to(device) #feature position in world frame\n",
    "        \n",
    "        self._V_flatten=torch.tensor([100.,1.,1.]).to(device)\n",
    "        self._V=torch.diag(self._V_flatten).to(device)\n",
    "        self._V_inv=torch.diag(1/self._V_flatten).to(device)\n",
    "        \n",
    "        #self._R_W_C=zy_rotation_matrix(angle_z=self._s[0],angle_y=self._s[1]) # camera pose in World frame\n",
    "        self._input_dim=max_num_landmarks*4+4+4+4\n",
    "        self._fov=fov\n",
    "        \n",
    "        self._policy=PolicyNetAtt(input_dim=self._input_dim).to(device)\n",
    "        self._policy_optimizer=Adam(self._policy.parameters(), lr=lr)\n",
    "        self.loss=0.\n",
    "    \n",
    "    def plan(self,s,p_b,p_T,p_f):\n",
    "        device=self._device\n",
    "        padding=torch.zeros(3*(self._max_num_landmarks-p_f.shape[1])).to(device)\n",
    "\n",
    "        \n",
    "        mask=torch.tensor([True]*(p_f.shape[1]+2)+[False]*(self._max_num_landmarks-p_f.shape[1])).to(device)\n",
    "        net_input=torch.hstack([s,p_b,p_T,p_f.flatten(),padding,mask]).to(device)\n",
    "\n",
    "        \n",
    "        action=self._policy.forward(net_input)\n",
    "        \n",
    "        return action \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "input dim:\n",
      "28\n",
      "torch.Size([1024, 3, 3])\n",
      "torch.Size([1024, 3])\n",
      "tensor([[ True,  True,  True,  True,  True, False],\n",
      "        [ True,  True,  True,  True,  True, False],\n",
      "        [ True,  True,  True,  True,  True, False],\n",
      "        ...,\n",
      "        [ True,  True,  True,  True,  True, False],\n",
      "        [ True,  True,  True,  True,  True, False],\n",
      "        [ True,  True,  True,  True,  True, False]], device='cuda:0')\n",
      "net_input.shape:\n",
      "torch.Size([1024, 28])\n",
      "estimated_landmark_pos.shape:\n",
      "torch.Size([1024, 18])\n",
      "mask.shape:\n",
      "torch.Size([1024, 1, 6])\n"
     ]
    }
   ],
   "source": [
    "max_landmark=4\n",
    "batch_size=1024\n",
    "\n",
    "ag=Agent(device=device,max_num_landmarks=max_landmark)\n",
    "ag._num_landmarks=3\n",
    "\n",
    "ag._s=torch.rand([batch_size,4]).to(device)\n",
    "ag._p_b=torch.rand([batch_size,3]).to(device)\n",
    "ag._p_T=torch.rand([batch_size,3]).to(device)\n",
    "\n",
    "\n",
    "\n",
    "ag._p_f=torch.rand([batch_size,3,ag._num_landmarks]).to(device)\n",
    "p_f_flatten=ag._p_f.reshape([batch_size,3*ag._num_landmarks]).to(device)\n",
    "\n",
    "padding=torch.zeros(3*(ag._max_num_landmarks-ag._p_f.shape[2])).to(device)\n",
    "padding=padding.expand(batch_size,padding.shape[0])\n",
    "\n",
    "mask=torch.tensor([True]*(ag._p_f.shape[2]+2)+[False]*(ag._max_num_landmarks-ag._p_f.shape[2])).to(device)\n",
    "mask=mask.expand(batch_size,mask.shape[0])\n",
    "\n",
    "\n",
    "print(ag._p_f.shape)\n",
    "print(padding.shape)\n",
    "print(mask)\n",
    "\n",
    "net_input=torch.hstack([ag._s,ag._p_b,ag._p_T,p_f_flatten,padding,mask])\n",
    "print(\"net_input.shape:\")\n",
    "print(net_input.shape)\n",
    "\n",
    "#net_input=torch.rand([batch_size,max_landmark*4+4+4+4]).to(device)\n",
    "action=ag._policy.forward(net_input)\n",
    "\n",
    "ag._s[:,:2]+=action\n",
    "\n",
    "R=euler_to_rotation_matrix_zyx(angle_z=ag._s[:,0],angle_y=ag._s[:,1],angle_x=ag._s[:,1]*0.)\n",
    "\n",
    "R=R.reshape([1024,3,3])\n",
    "\n",
    "ag._p_T.shape\n",
    "\n",
    "p_T_C=R@(ag._p_T.reshape([1024,3,1]))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "torch.Size([1024, 3, 1])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R = torch.rand(1024,3,3)\n",
    "v = torch.rand(1024,3,1)\n",
    "\n",
    "v_rotated=R@v\n",
    "dist=torch.sqrt(torch.sum(v_rotated*v_rotated,axis=1))\n",
    "v_x=v_rotated[:,0,:]\n",
    "bearing=v_x/dist\n",
    "\n",
    "v_rotated.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "torch.Size([1024, 3])\n"
     ]
    }
   ],
   "source": [
    "\n",
    "batch_size=1024\n",
    "yaw=torch.zeros([batch_size]).to(device)+torch.rand([batch_size]).to(device)*0.001\n",
    "pitch=torch.zeros([batch_size]).to(device)+torch.rand([batch_size]).to(device)*0.001\n",
    "\n",
    "R=euler_to_rotation_matrix_zyx(angle_z=yaw,angle_y=pitch,angle_x=yaw)\n",
    "R.shape\n",
    "\n",
    "\n",
    "R=torch.transpose(R,0,2)\n",
    "R=torch.transpose(R,1,2)\n",
    "\n",
    "\n",
    "pts = torch.rand(1024,3,6).to(device)\n",
    "\n",
    "pts_rotated=R@pts\n",
    "\n",
    "i=2\n",
    "    \n",
    "v=pts_rotated[:,:,i]*0+torch.tensor([1.,-0.1,0.]).to(device)\n",
    "v_norm=torch.sqrt(torch.sum(v*v,axis=1))\n",
    "v=v/(v_norm.repeat(3,1).T)\n",
    "\n",
    "\n",
    "\n",
    "unit_x=torch.tensor([1.,0.,0.]).to(device).repeat(batch_size,1)\n",
    "\n",
    "\n",
    "c=torch.cross(unit_x,v)\n",
    "\n",
    "dot=torch.sum(unit_x * v,axis=1)\n",
    "\n",
    "V=torch.zeros([1024,3,3]).to(device)\n",
    "V[:,0,1]=-c[:,2]\n",
    "V[:,0,2]=c[:,1]\n",
    "V[:,1,0]=c[:,2]\n",
    "V[:,1,2]=-c[:,0]\n",
    "V[:,2,0]=-c[:,1]\n",
    "V[:,2,1]=c[:,0]\n",
    "\n",
    "R_W_C=torch.eye(3).repeat(1024,1,1).to(device)\n",
    "\n",
    "\n",
    "c_norm=torch.sqrt(torch.sum(c*c,axis=1))\n",
    "\n",
    "weight=((1-dot)/(c_norm**2)).repeat(3,3,1)\n",
    "weight=torch.transpose(weight,0,2).to(device)\n",
    "\n",
    "R_W_C=R_W_C+V+weight*V@V\n",
    "R_C_W=torch.transpose(R_W_C,1,2)\n",
    "\n",
    "\n",
    "\n",
    "V_inv=torch.diag(torch.tensor([40.,1.,1.])).to(device)\n",
    "R_W_C.shape\n",
    "\n",
    "V_inv=V_inv.repeat(1024,1,1)\n",
    "\n",
    "\n",
    "\n",
    "S=R_W_C @ V_inv @ R_C_W"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "最优解： 1.0\n"
     ]
    }
   ],
   "source": [
    "import torch\n",
    "import torch.optim as optim\n",
    "\n",
    "# 定义一个待优化的变量\n",
    "x = torch.tensor([1.0], requires_grad=True)\n",
    "\n",
    "# 定义损失函数，这里以简单的二次函数为例\n",
    "def loss_function(x):\n",
    "    return 2*x**2 - 4*x + 1\n",
    "\n",
    "# 创建一个优化器，这里使用梯度下降法\n",
    "optimizer = optim.SGD([x], lr=0.1)\n",
    "\n",
    "# 进行优化迭代\n",
    "for i in range(100):\n",
    "    optimizer.zero_grad()  # 梯度清零\n",
    "    loss = loss_function(x)\n",
    "    loss.backward()  # 反向传播计算梯度\n",
    "    optimizer.step()  # 更新变量\n",
    "\n",
    "print(\"最优解：\", x.item())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor(4.5763, grad_fn=<MulBackward0>)\n",
      "tensor(9.2173, grad_fn=<MulBackward0>)\n",
      "tensor(9.8731, grad_fn=<MulBackward0>)\n",
      "tensor(9.9714, grad_fn=<MulBackward0>)\n",
      "tensor(9.9919, grad_fn=<MulBackward0>)\n",
      "tensor(9.9971, grad_fn=<MulBackward0>)\n",
      "tensor(9.9988, grad_fn=<MulBackward0>)\n",
      "tensor(9.9995, grad_fn=<MulBackward0>)\n",
      "tensor(9.9997, grad_fn=<MulBackward0>)\n",
      "tensor(9.9999, grad_fn=<MulBackward0>)\n",
      "tensor(9.9999, grad_fn=<MulBackward0>)\n",
      "tensor(10.0000, grad_fn=<MulBackward0>)\n",
      "tensor(10.0000, grad_fn=<MulBackward0>)\n",
      "tensor(10., grad_fn=<MulBackward0>)\n",
      "tensor(10.0000, grad_fn=<MulBackward0>)\n",
      "tensor(10.0000, grad_fn=<MulBackward0>)\n",
      "tensor(10.0000, grad_fn=<MulBackward0>)\n",
      "tensor(10.0000, grad_fn=<MulBackward0>)\n",
      "tensor(10.0000, grad_fn=<MulBackward0>)\n",
      "tensor(9.9999, grad_fn=<MulBackward0>)\n",
      "0.24664831161499023\n",
      "tensor([-0.6336,  0.6563], requires_grad=True)\n",
      "tensor(-110.4639, grad_fn=<SubBackward0>)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f38f166eb80>]"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjMAAAGdCAYAAADnrPLBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABBJElEQVR4nO3deXxU9aH+8WeWzGQhG2QDCZAAAkGkLIpBsFQRVLB622qtFcS9XtSqVIGfCi61aG1p1dqqvbLYel26WL24YFBQlgjKTtgFErYsbJmEkExm5vz+CBkTSYBAJmeWz/v1mhfMOd85eY6nQ56eOXO+FsMwDAEAAIQoq9kBAAAAzgZlBgAAhDTKDAAACGmUGQAAENIoMwAAIKRRZgAAQEijzAAAgJBGmQEAACHNbnaAtuDz+bRv3z7Fx8fLYrGYHQcAAJwGwzBUUVGhTp06yWpt/vxLRJSZffv2KTMz0+wYAADgDOzevVudO3dudn1ElJn4+HhJdf8xEhISTE4DAABOh8vlUmZmpv/3eHMioszUf7SUkJBAmQEAIMSc6hIRLgAGAAAhjTIDAABCGmUGAACENMoMAAAIaZQZAAAQ0igzAAAgpFFmAABASKPMAACAkEaZAQAAIY0yAwAAQhplBgAAhDTKDAAACGmUmbPw0fr9evDtNSrYV252FAAAIpapZaZbt26yWCyNHs8880yjMevWrdPw4cMVHR2tzMxM/fa3vzUp7YneXb1X/169V4u2lJkdBQCAiGX6mZknn3xS+/fv9z/uvfde/zqXy6VRo0apa9euWrlypZ577jk9/vjjevXVV01M/K3hPVMkSYu3UWYAADCL3ewA8fHxysjIaHLdG2+8IbfbrVmzZsnhcKhv375as2aNZs6cqTvvvLONk55oWM9USdLKwsOqcnsU6zD9PycAABHH9DMzzzzzjDp06KABAwboueeek8fj8a/Lz8/XJZdcIofD4V82evRobdmyRYcPHzYjbiPdOsTqnKQY1XoNLd95yOw4AABEJFNPJdx3330aOHCg2rdvr2XLlmnq1Knav3+/Zs6cKUkqLi5WVlZWo9ekp6f71yUnJze53ZqaGtXU1Pifu1yugOS3WCwa3jNFb321W0u2HdAPeqUF5OcAAIDmtfqZmSlTppxwUe93H5s3b5YkPfjggxoxYoTOP/98/eIXv9Dvf/97vfjii42KyJmYMWOGEhMT/Y/MzMzW2LUmDTt+3cySbQcC9jMAAEDzWv3MzKRJkzRhwoSTjsnOzm5y+ZAhQ+TxeLRr1y716tVLGRkZKikpaTSm/nlz19lI0tSpU/Xggw/6n7tcroAVmqHdU2SxSFtKKlTqqlZaQnRAfg4AAGhaq5eZ1NRUpaamntFr16xZI6vVqrS0uo9rcnNz9cgjj6i2tlZRUVGSpLy8PPXq1avZj5gkyel0yul0nlGGlmof51DfTgnasNelJdsP6EcDO7fJzwUAAHVMuwA4Pz9ff/zjH7V27Vrt2LFDb7zxhh544AHddNNN/qJy4403yuFw6LbbblNBQYHefvttPf/8843OugSDYT3qytuS7XzUBABAWzPtAmCn06m33npLjz/+uGpqapSVlaUHHnigUVFJTEzUJ598ookTJ2rQoEFKSUnRtGnTguJr2Q0N75milz//Rku2HZBhGLJYLGZHAgAgYphWZgYOHKgvv/zylOPOP/98LV68uA0SnblBXZPltFtVWlGjbaWVOjc93uxIAABEDNPvMxMOoqNsujCrvSRpMd9qAgCgTVFmWslw/1e0mdoAAIC2RJlpJfUXAS/feUhuj8/kNAAARA7KTCvpnRGvlHYOVbm9WlVk/lQLAABECspMK7FaLRranbsBAwDQ1igzrah+aoPF3G8GAIA2Q5lpRfUXAa/fc0TlVbUmpwEAIDJQZlpRx8QYdU+Nk8+Q8ndwdgYAgLZAmWllw3vWfauJ+80AANA2KDOtbFiP4xcBc90MAABtgjLTyi7q3kF2q0WFB6u0+1CV2XEAAAh7lJlW1s5p14AuSZL4qAkAgLZAmQmA+rsBL9nO1AYAAAQaZSYAhvXsIElauv2gvD7D5DQAAIQ3ykwA9O+cpHinXeXHarVhb7nZcQAACGuUmQCw26y6qHvd2Rm+1QQAQGBRZgKk/m7AzNMEAEBgUWYCpP5+MysLD+uY22tyGgAAwhdlJkCyUuJ0TlKM3F6flu88aHYcAADCFmUmQCwWy7d3A+ajJgAAAoYyE0DDejK1AQAAgUaZCaCLj5+Z2VxcodKKapPTAAAQnigzAdQ+zqG+nRIkSUs5OwMAQEBQZgKs/qMm5mkCACAwKDMBNrx+nqZtB2QYTG0AAEBro8wE2OBuyXLarSqtqNG20kqz4wAAEHYoMwEWHWXThVntJfEVbQAAAoEy0wb895vhImAAAFodZaYN1F8E/OWOg3J7fCanAQAgvFBm2kCfjAR1iHOoyu3V6qLDZscBACCsUGbagNVq0VA+agIAICAoM21keA/uNwMAQCBQZtpI/XUz6/YcUXlVrclpAAAIH5SZNtIpKUbZqXHyGVL+Ds7OAADQWigzbYiPmgAAaH2UmTY0rOfxqQ24CBgAgFZDmWlDF2W3l81qUeHBKu0+VGV2HAAAwgJlpg3FR0dpQGaSJM7OAADQWigzbaz+W03M0wQAQOugzLSx+nmaln5zQF6fYXIaAABCH2WmjfXPTFI7p11HqmpVsK/c7DgAAIQ8ykwbi7JZdVF2B0l8RRsAgNZAmTHBcK6bAQCg1VBmTFB/EfDKwsM65vaanAYAgNBmepn54IMPNGTIEMXExCg5OVnXXntto/VFRUUaM2aMYmNjlZaWpoceekgej8ecsK0kOyVOnRKj5fb6tHznQbPjAAAQ0uxm/vB//etfuuOOO/Sb3/xGl156qTwejzZs2OBf7/V6NWbMGGVkZGjZsmXav3+/xo8fr6ioKP3mN78xMfnZsVgsGtYzRe98vUdLtx/QiF5pZkcCACBkWQzDMOX7wR6PR926ddMTTzyh2267rckxH330kcaOHat9+/YpPT1dkvTyyy9r8uTJKisrk8PhOK2f5XK5lJiYqPLyciUkJLTaPpyN99fu031vrlbvjHh9fP8lZscBACDonO7vb9M+Zlq1apX27t0rq9WqAQMGqGPHjrryyisbnZnJz89Xv379/EVGkkaPHi2Xy6WCggIzYreaod3rvtG0ubhCZRU1JqcBACB0mVZmduzYIUl6/PHH9eijj2revHlKTk7WiBEjdOjQIUlScXFxoyIjyf+8uLi42W3X1NTI5XI1egSblHZO5XSsa5lLmdoAAIAz1uplZsqUKbJYLCd9bN68WT6fT5L0yCOP6Mc//rEGDRqk2bNny2Kx6B//+MdZZZgxY4YSExP9j8zMzNbYtVZX/xVt7jcDAMCZa/ULgCdNmqQJEyacdEx2drb2798vScrJyfEvdzqdys7OVlFRkSQpIyNDK1asaPTakpIS/7rmTJ06VQ8++KD/ucvlCspCM6xnil75YoeWbC+TYRiyWCxmRwIAIOS0eplJTU1VamrqKccNGjRITqdTW7Zs0bBhwyRJtbW12rVrl7p27SpJys3N1dNPP63S0lKlpdV94ycvL08JCQmNStB3OZ1OOZ3OVtibwLqgW3s57FaVuGq0vbRSPdPjzY4EAEDIMe2amYSEBP3iF7/Q9OnT9cknn2jLli26++67JUnXXXedJGnUqFHKycnRuHHjtHbtWs2fP1+PPvqoJk6cGBJl5VSio2y6sFt7SXzUBADAmTL1pnnPPfecbrjhBo0bN04XXHCBCgsL9dlnnyk5OVmSZLPZNG/ePNlsNuXm5uqmm27S+PHj9eSTT5oZu1XV3w14CRcBAwBwRky7z0xbCsb7zNTbsLdcY19coliHTWumjZLDbvpNmQEACApBf58Z1MnpmKD2cQ5Vub1as/uI2XEAAAg5lBmTWa0W/w30lmwrMzkNAAChhzITBPz3m+G6GQAAWowyEwSG9az7Kvva3UdUfqzW5DQAAIQWykwQOCcpRtkpcfIZUv43B82OAwBASKHMBIlvv6LNdTMAALQEZSZIDOtxvMxw8zwAAFqEMhMkLureQTarRbsOVmn3oSqz4wAAEDIoM0EiITpK38tMksTdgAEAaAnKTBC5uAdTGwAA0FKUmSBSf7+ZZdsPyOcL+1kmAABoFZSZIPK9zCS1c9p1uKpWBftcZscBACAkUGaCSJTNqouy20uSFvMVbQAATgtlJsjwFW0AAFqGMhNk6qc2+HrXYR1ze01OAwBA8KPMBJnuqXHqmBgtt9enFbsOmR0HAICgR5kJMhaLpcFHTVw3AwDAqVBmglD9PE2LuW4GAIBToswEofqb520urlBZRY3JaQAACG6UmSCU0s6pPh0TJEnLvuHsDAAAJ0OZCVLD+agJAIDTQpkJUg3vN2MYTG0AAEBzKDNB6sKs9nLYrSp2VWvdnnKz4wAAELQoM0EqOsqmq87LkCT9dfEOk9MAABC8KDNB7M5LukuSPly/X0UHq0xOAwBAcKLMBLGcTgm65NxU+Qzpf5ZwdgYAgKZQZoLcLy7JliS98/VuHazknjMAAHwXZSbI5XbvoH7nJKq61qfX8wvNjgMAQNChzAQ5i8Wiu75fd3bm9fxdqnJ7TE4EAEBwocyEgCv6ZqhL+1gdrqrVP77eY3YcAACCCmUmBNhtVt0xPEtS3de0PV6fyYkAAAgelJkQ8ZNBmWof59Cew8f04YZis+MAABA0KDMhIsZh08253SRJr3z+DVMcAABwHGUmhIzP7aqYKJsK9rm0dPtBs+MAABAUKDMhJDnOoZ9ekClJeuWLb0xOAwBAcKDMhJjbhmXJZrVo8bYD2rCXCSgBAKDMhJjM9rEa06+jJOnVL5jiAAAAykwIuvP4FAcfrN+v3YeYgBIAENkoMyHovHMSNbxnirw+Q68t2Wl2HAAATEWZCVF3XdJdkvT2V7t1+Kjb5DQAAJiHMhOiLu7RQX07JehYrVd/+5IJKAEAkYsyE6LqJqCsOzszZ9kuVdd6TU4EAIA5KDMh7KrzMtQ5OUaHjrr1j5VMQAkAiEyUmRBWNwFl3Teb/vrFDnl9THEAAIg8ppWZRYsWyWKxNPn46quv/OPWrVun4cOHKzo6WpmZmfrtb39rVuSgdN3gzkqOjVLRoSp9zASUAIAIZFqZGTp0qPbv39/ocfvttysrK0uDBw+WJLlcLo0aNUpdu3bVypUr9dxzz+nxxx/Xq6++albsoBPrsGv88QkoX2YCSgBABDKtzDgcDmVkZPgfHTp00HvvvadbbrlFFotFkvTGG2/I7XZr1qxZ6tu3r2644Qbdd999mjlzplmxg9L43K6KjrJq/d5y5e9gAkoAQGQJmmtm3n//fR08eFC33HKLf1l+fr4uueQSORwO/7LRo0dry5YtOnz4cLPbqqmpkcvlavQIZx3aOXX94OMTUH7OFAcAgMgSNGXmtdde0+jRo9W5c2f/suLiYqWnpzcaV/+8uLj560NmzJihxMRE/yMzMzMwoYPI7cOyZbVIn28t06b94V3eAABoqNXLzJQpU5q9sLf+sXnz5kav2bNnj+bPn6/bbrutVTJMnTpV5eXl/sfu3btbZbvBrEuHWF3FBJQAgAhkb+0NTpo0SRMmTDjpmOzs7EbPZ8+erQ4dOuiHP/xho+UZGRkqKSlptKz+eUZGRrPbdzqdcjqdLUgdHu66pLvmrduv99fu06RR56pzcqzZkQAACLhWLzOpqalKTU097fGGYWj27NkaP368oqKiGq3Lzc3VI488otraWv+6vLw89erVS8nJya2aOxz065yoi3t00NLtBzVryS5NuzrH7EgAAASc6dfMfPbZZ9q5c6duv/32E9bdeOONcjgcuu2221RQUKC3335bzz//vB588EETkoaG+gko3/qqSEeqmIASABD+TC8zr732moYOHarevXufsC4xMVGffPKJdu7cqUGDBmnSpEmaNm2a7rzzThOShobhPVPUp2OCqtxe/Z0JKAEAEcBiRMBd1lwulxITE1VeXq6EhASz4wTce2v26pdvrVFKO4eWTL5U0VE2syMBANBip/v72/QzM2h9V/XrqHOSYnSg0q1/rWICSgBAeKPMhKEom1W3D8+SxASUAIDwR5kJU9cPzlRiTJR2HazSJwVMQAkACF+UmTAV57RrfG5XSUxACQAIb5SZMHbz0G5y2K1au6dcy3ceMjsOAAABQZkJYyntnLpuUN1cV698/o3JaQAACAzKTJi7Y3i2LBZp4ZYybS5mAkoAQPihzIS5bilxuvK8unmsmIASABCOKDMRoH6Kg/fX7NO+I8dMTgMAQOuizESA/plJuii7vTw+Q7OW7DQ7DgAArYoyEyHu+n7d2Zk3VxSpvKrW5DQAALQeykyEGHFuqnqlx+uo26u/L2cCSgBA+KDMRAiLxaK7vp8tSZq9dJeqa70mJwIAoHVQZiLI1f07qVNitA5U1ug/q/eaHQcAgFZBmYkgUTarbh1WNwHlq1/skI8JKAEAYYAyE2FuuLCLEqLt2nHgqPI2lZgdBwCAs0aZiTDtnHaNYwJKAEAYocxEoPoJKFcXHdHXhYfNjgMAwFmhzESgtPho/Xhg3QSUf1nEBJQAgNBGmYlQd16SLZvVos82l2rp9gNmxwEA4IxRZiJUVkqcxl1Ud+3ME/9XoFqvz+REAACcGcpMBHtg5LlKjo3S1pJK/f1L7goMAAhNlJkIlhgbpV+N7iVJ+kPeVh2srDE5EQAALUeZiXA3XNBFOR0T5Kr26HefbDU7DgAALUaZiXA2q0VPXNNXkvTWV0XasLfc5EQAALQMZQa6oFt7/bB/JxmG9Pj7BdxIDwAQUigzkCRNvaq3YqJs+rrwsN5fu8/sOAAAnDbKDCRJHRNjdM+lPSRJv/lwk47WeExOBADA6aHMwO+2YVnq0j5WJa4a/XnRdrPjAABwWigz8IuOsunRMX0kSX/9YqcKDx41OREAAKdGmUEjl+eka3jPFLm9Pj01b5PZcQAAOCXKDBqxWCyafnWO7FaLFmwq0edby8yOBADASVFmcIIeafG6eWg3SXXzNrk9zNsEAAhelBk06ZcjeyqlnUM7yo7q9fxdZscBAKBZlBk0KSE6Sg+P7i1Jen7BNpVVMG8TACA4UWbQrJ8M6qzzOyeqosaj5+ZvNjsOAABNosygWVarRdOvrpu36Z2v92jN7iPmBgIAoAmUGZzUoK7J+tHAcyTVzdvk8zFvEwAguFBmcEpTruitOIdNa3Yf0bur95odBwCARigzOKW0hGjde1lPSdIzH29WRXWtyYkAAPgWZQan5ZaLuykrJU5lFTX602fM2wQACB6UGZwWp92mx8bWzds0a+lO7SirNDkRAAB1KDM4bZf2TtcPeqWq1mvoqXkbzY4DAIAkygxa6LGxOYqyWbRwS5k+21xidhwAAMwtM1u3btU111yjlJQUJSQkaNiwYVq4cGGjMUVFRRozZoxiY2OVlpamhx56SB6Px6TEyE5tp1uHZUmSnpq3STUer8mJAACRztQyM3bsWHk8Hn322WdauXKl+vfvr7Fjx6q4uFiS5PV6NWbMGLndbi1btkxz587VnDlzNG3aNDNjR7x7L+2p1Hindh44qtlLd5kdBwAQ4SyGYZhyF7QDBw4oNTVVX3zxhYYPHy5JqqioUEJCgvLy8jRy5Eh99NFHGjt2rPbt26f09HRJ0ssvv6zJkyerrKxMDofjtH6Wy+VSYmKiysvLlZCQELB9iiT/WrlHk/6xVnEOmxb+aoTSEqLNjgQACDOn+/vbtDMzHTp0UK9evfT666/r6NGj8ng8euWVV5SWlqZBgwZJkvLz89WvXz9/kZGk0aNHy+VyqaCgoNlt19TUyOVyNXqgdf3XgHM0oEuSjrq9euZj5m0CAJjHtDJjsVi0YMECrV69WvHx8YqOjtbMmTP18ccfKzk5WZJUXFzcqMhI8j+v/yiqKTNmzFBiYqL/kZmZGbgdiVBWq0WPH5+36d+r9mpl4WGTEwEAIlWrl5kpU6bIYrGc9LF582YZhqGJEycqLS1Nixcv1ooVK3Tttdfq6quv1v79+88qw9SpU1VeXu5/7N69u5X2Dg31z0zS9YM7S5Ke+D/mbQIAmMPe2hucNGmSJkyYcNIx2dnZ+uyzzzRv3jwdPnzY/znYn//8Z+Xl5Wnu3LmaMmWKMjIytGLFikavLSmp+zpwRkZGs9t3Op1yOp1ntyM4LQ+N7q2P1hdr3Z5y/XPlHl1/AWfBAABtq9XLTGpqqlJTU085rqqqSpJktTY+OWS1WuXz+SRJubm5evrpp1VaWqq0tDRJUl5enhISEpSTk9PKyXEmUuOd+uXInvr1B5v07MebNfq8DCXGRJkdCwAQQUy7ZiY3N1fJycm6+eabtXbtWm3dulUPPfSQdu7cqTFjxkiSRo0apZycHI0bN05r167V/Pnz9eijj2rixImceQki43O7qXtqnA4edeuFT7eZHQcAEGFMKzMpKSn6+OOPVVlZqUsvvVSDBw/WkiVL9N5776l///6SJJvNpnnz5slmsyk3N1c33XSTxo8fryeffNKs2GiCw27V9OMXA89dtkvbSytMTgQAiCSm3WemLXGfmbZx+9yvtWBTiYb3TNHrt14oi8VidiQAQAgL+vvMIPw8NraPHDarFm87oLyNzNsEAGgblBm0mq4d4nTHJXXzNv36g02qrmXeJgBA4FFm0Kr+e0QPZSREq+hQlV5bstPsOACACECZQauKc9o19arekqQ/fbZd+8uPmZwIABDuKDNodT/s30mDuybrWK1Xj79foAi4xhwAYCLKDFqdxWLRk9ecJ7vVovkFJXp/7T6zIwEAwhhlBgGR0ylBv7yspyTpsf9sUImr2uREAIBwRZlBwNw9orv6d06Uq9qjyf9ax8dNAICAoMwgYOw2q35/fX857FYt2lKmt79i9nIAQOujzCCgeqTF6+HRvSRJT83bqN2HqkxOBAAIN5QZBNwtF2fpgm7JOur26uF/rpPPx8dNAIDWQ5lBwNmsFv3uuv6KibIpf8dBvZ6/y+xIAIAwQplBm+jaIU7/b0wfSdIzH2/WjrJKkxMBAMIFZQZt5qYhXTS8Z4qqa32a9I+18vJxEwCgFVBm0GYsFoue/fH5infatbroiF79YofZkQAAYYAygzbVKSlG067OkST9IW+rthRXmJwIABDqKDNocz8Z1Fkj+6TJ7fXpwXfWqNbrMzsSACCEUWbQ5iwWi37zo35Kio1SwT6X/vTZdrMjAQBCGGUGpkiLj9avrz1PkvSnhdu1fk+5yYkAAKGKMgPTjD2/k8ae31Fen6EH31mj6lqv2ZEAACGIMgNTPXXNeUpp59S20kr9YcFWs+MAAEIQZQamSo5zaMaP+kmSXv1ih1YWHjI5EQAg1FBmYLrLc9L1k0GdZRjSpHfWqsrtMTsSACCEUGYQFKZdnaNOidHadbBKz3602ew4AIAQQplBUEiIjtJvf9JfkjQ3v1BLtx8wOREAIFRQZhA0hvVM0U0XdZEkPfzPdaqorjU5EQAgFFBmEFSmXtlHXdrHau+RY/r1vE1mxwEAhADKDIJKnNOu313XXxaL9PbXu/XZ5hKzIwEAghxlBkHnwqz2un1YliRp8r/W6/BRt8mJAADBjDKDoDRpVC91T41TWUWNpr9fYHYcAEAQo8wgKEVH2fT7678nm9Wi99fu04fr95sdCQAQpCgzCFrfy0zSf4/oLkl69D8bVFZRY3IiAEAwoswgqN17aU/16ZigQ0fdeuTd9TIMw+xIAIAgQ5lBUHPYrZp5fX9F2Sz6ZGOJ3l291+xIAIAgQ5lB0OvTMUH3jzxXkjT9/QLtLz9mciIAQDChzCAk3HVJtvpnJqmi2qPJ/+LjJgDAtygzCAl2m1W/v66/nHarvthapjdX7DY7EgAgSFBmEDJ6pLXTw1f0liT9+oONKjpYZXIiAEAwoMwgpNwytJuGZLVXldurX/1zrXw+Pm4CgEhHmUFIsVoteu4n/RXrsGnFzkN64bNtZkcCAJiMMoOQ06VDrJ685jxJ0h8XbNP8gmKTEwEAzESZQUj6yaDOmjC0myTpwbfXaEtxhbmBAACmocwgZD0ypo+Gdu+go26v7nj9ax2pYnZtAIhElBmErCibVS/dOFCZ7WNUdKhK9/zvanm8PrNjAQDamKllZtWqVbr88suVlJSkDh066M4771RlZWWjMUVFRRozZoxiY2OVlpamhx56SB6Px6TECDbJcQ79dfxgxTpsWrL9gGZ8tNnsSACANmZamdm3b59GjhypHj16aPny5fr4449VUFCgCRMm+Md4vV6NGTNGbrdby5Yt09y5czVnzhxNmzbNrNgIQr0zEjTz+v6SpNeW7NQ/V+4xOREAoC1ZDJPuC//qq6/qscce0/79+2W11nWq9evX6/zzz9e2bdvUo0cPffTRRxo7dqz27dun9PR0SdLLL7+syZMnq6ysTA6H47R+lsvlUmJiosrLy5WQkBCwfYK5/pC3Vc9/uk0Ou1Vv33mRBnRJNjsSAOAsnO7vb9POzNTU1MjhcPiLjCTFxMRIkpYsWSJJys/PV79+/fxFRpJGjx4tl8ulgoKCk27b5XI1eiD8/fKynhqVky63x6e7/rZSJa5qsyMBANqAaWXm0ksvVXFxsZ577jm53W4dPnxYU6ZMkSTt379fklRcXNyoyEjyPy8ubv7eIjNmzFBiYqL/kZmZGaC9QDCxWi2a+dPv6dz0diqtqNFdf1up6lqv2bEAAAHW6mVmypQpslgsJ31s3rxZffv21dy5c/X73/9esbGxysjIUFZWltLT0xudrTkTU6dOVXl5uf+xezeTEkaKdk67/jp+sBJjorRm9xE98u4GZtgGgDBnb+0NTpo0qdFFvE3Jzs6WJN1444268cYbVVJSori4OFksFs2cOdO/PiMjQytWrGj02pKSEv+65jidTjmdzrPYC4Syrh3i9NKNAzV+1nL9a9Ue9e2UoFuHZZkdCwAQIK1eZlJTU5Wamtqi19R/dDRr1ixFR0fr8ssvlyTl5ubq6aefVmlpqdLS0iRJeXl5SkhIUE5OTusGR1gZ1jNFj4zJ0VPzNurpDzfp3PR4DeuZYnYsAEAAmHqfmT/96U9atWqVtm7dqpdeekn33HOPZsyYoaSkJEnSqFGjlJOTo3Hjxmnt2rWaP3++Hn30UU2cOJEzLzilWy/uph8P7Cyvz9DE/12lwoNHzY4EAAgAU8vMihUrdPnll6tfv3569dVX9corr+i+++7zr7fZbJo3b55sNptyc3N10003afz48XryySdNTI1QYbFY9PR/nafvZSap/Fit7nj9a1XWcMNFAAg3pt1npi1xn5nIVuKq1tUvLlFpRY1G903XX34+SFarxexYAIBTCPr7zABtJT0hWq+MGySHzar5BSV64bNtZkcCALQiygwiwoAuyfrNj/pJkv64YJs+3tD8fYoAAKGFMoOI8ZNBnXXrxXVf0X7wnTXaXMydoQEgHFBmEFH+31W9dXGPDqpye3XH61/r8FG32ZEAAGeJMoOIYrdZ9aefDVSX9rHafeiY7nlzlTxen9mxAABngTKDiJMc59Bfxw9WrMOmpdsP6ukPN5kdCQBwFigziEi9MuI18/rvSZJmL92ld75m/i4ACFWUGUSsK87L0AMjz5UkPfruBq0qOmxyIgDAmaDMIKLde2kPXdE3Q26vT7/420qVuKrNjgQAaCHKDCKa1WrR76/vr17p8SqtqNGdf1up6lqv2bEAAC1AmUHEi3Pa9dfxg5UUG6W1u4/okXc3KAJm+QCAsEGZASR16RCrl24cKJvVon+t2qNZS3eZHQkAcJooM8BxF/dI0aNj+kiSnv5go5ZsO2ByIgDA6aDMAA1MGNpN1w3qLJ8h/fcbK7V+T7nZkQAAp0CZARqwWCz69X+dpwu6JctV7dGN//Ol1u4+YnYsAMBJUGaA73DabZp9y4W6oFuyKqo9uul/lnMPGgAIYpQZoAntnHbNueVCXZjVXhU1Ho1/bYVWFh4yOxYAoAmUGaAZcU675txygXKzO6jyeKFZsZNCAwDBhjIDnESsw65ZEy7QsB4pOur26uZZK5T/zUGzYwEAGqDMAKcQ47Dpf24erOE9U3Ss1qtb5qzQ0u18bRsAggVlBjgN0VE2/XX8YI3olarqWp9unfOVvthaZnYsAIAoM8Bpi46y6ZVxg3RZ7zTVeHy6/fWvtWhLqdmxACDiUWaAFnDabfrLTYN0eU663B6f7nx9pRZuptAAgJkoM0ALOexW/fnnA3XleRlye326829fa8HGErNjAUDEoswAZyDKZtULPxugMf06qtZr6O43VurjDcVmxwKAiESZAc5QlM2q52/4nq7u30m1XkP3/O8qfbh+v9mxACDiUGaAs2C3WfWH6/vrvwacI4/P0L1vrtb/rd1ndiwAiCiUGeAs2W1W/e66/vrxwM7y+gz98q3Vem/NXrNjAUDEoMwArcBmtei5n5yvnw7OlM+QHnh7jf61co/ZsQAgIlBmgFZitVo040f99LMLu8hnSL/651q98/Vus2MBQNijzACtyGq16Olrz9NNF3WRYUiT/7VOb60oMjsWAIQ1ygzQyqxWi5665jxNGNpNhiFN+fd6vbG80OxYABC2KDNAAFgsFk2/Oke3XpwlSXrk3Q16PX+XuaEAIExRZoAAsVgsemxsH915SbYkadp7BZq1ZKfJqQAg/FBmgACyWCyaemVv3T2iuyTpyXkb9T+Ld5icCgDCC2UGCDCLxaKHR/fSvZf2kCT9+oNNevnzb0xOBQDhgzIDtAGLxaJJo3rp/pE9JUnPfLRZMz/ZIp/PMDkZAIQ+ygzQhu4fea4mXX6uJOmFz7brljlf6dBRt8mpACC0UWaANnbvZT3125+cL6fdqs+3lmnsC4u1quiw2bEAIGRRZgATXD84U/+ZeLGyUuK0r7xa17+cr1lLdsow+NgJAFqKMgOYpE/HBL1/z8Ua06+jPD5DT87bqP9+Y5Vc1bVmRwOAkEKZAUwUHx2lP904QI9fnaMom0UfbSjWD19cooJ95WZHA4CQQZkBTGaxWDTh4iy9c1euzkmK0a6DVfqvPy/TWyuK+NgJAE4DZQYIEgO6JGvevcP0g16pcnt8mvLv9Zr0j7WqcnvMjgYAQS2gZebpp5/W0KFDFRsbq6SkpCbHFBUVacyYMYqNjVVaWpoeeugheTyN//FetGiRBg4cKKfTqR49emjOnDmBjA2YJjnOodduvkAPje4lq0X696q9uvalpdpeWml2NAAIWgEtM263W9ddd53uvvvuJtd7vV6NGTNGbrdby5Yt09y5czVnzhxNmzbNP2bnzp0aM2aMfvCDH2jNmjW6//77dfvtt2v+/PmBjA6Yxmq1aOIPeuh/77hIqfFObS2p1DV/WqL31+4zOxoABCWL0QYfys+ZM0f333+/jhw50mj5Rx99pLFjx2rfvn1KT0+XJL388suaPHmyysrK5HA4NHnyZH3wwQfasGGD/3U33HCDjhw5oo8//vi0fr7L5VJiYqLKy8uVkJDQavsFBFppRbXue3O1vtxxSJI07qKuenRsHzntNpOTAUDgne7vb1OvmcnPz1e/fv38RUaSRo8eLZfLpYKCAv+YkSNHNnrd6NGjlZ+f3+x2a2pq5HK5Gj2AUJQWH62/3zZEE39QN1Hl374s1HUv52v3oSqTkwFA8DC1zBQXFzcqMpL8z4uLi086xuVy6dixY01ud8aMGUpMTPQ/MjMzA5AeaBt2m1UPje6t2RMuUFJslNbtKdeYFxZrwcYSs6MBQFBocZmZMmWKLBbLSR+bN28ORNbTNnXqVJWXl/sfu3fvNjUP0Bp+0DtNH9w3XN/LTJKr2qPbX/9aMz7aJI/XZ3Y0ADCVvaUvmDRpkiZMmHDSMdnZ2ae1rYyMDK1YsaLRspKSEv+6+j/rlzUck5CQoJiYmCa363Q65XQ6TysDEErOSYrRO3flasZHmzR76S698vkOrS48ohdvHKD0hGiz4wGAKVpcZlJTU5WamtoqPzw3N1dPP/20SktLlZaWJknKy8tTQkKCcnJy/GM+/PDDRq/Ly8tTbm5uq2QAQo3DbtX0q/vqgm7t9fA/12nFrkMa88JiPX/DAF3cI8XseADQ5gJ6zUxRUZHWrFmjoqIieb1erVmzRmvWrFFlZd09M0aNGqWcnByNGzdOa9eu1fz58/Xoo49q4sSJ/jMrv/jFL7Rjxw49/PDD2rx5s/785z/rnXfe0QMPPBDI6EDQu6pfR/3fvcPUOyNeByrduum15Xrh023y+bhrMIDIEtCvZk+YMEFz5849YfnChQs1YsQISVJhYaHuvvtuLVq0SHFxcbr55pv1zDPPyG7/9qTRokWL9MADD2jjxo3q3LmzHnvssVN+1NUQX81GOKuu9Wr6ewV6++u6a8MuOTdVf/zp99Q+zmFyMgA4O6f7+7tN7jNjNsoMIsE/V+7Ro/9Zr+panzomRmvGj/ppRK80s2MBwBkLifvMAGg9PxnUWf+ZeLGyU+K0v7xaE2Z/pVtmr2AqBABhjzIDhJHeGQl6/95hun1YluxWixZuKdMVf/xCT/xfgcqras2OBwABwcdMQJjaUVap33y4SQs2lUqSkmOj9ODl5+pnF3aR3cb/jwEQ/LhmpgHKDCLZ4m1lemreRm0tqfu46dz0dnpsbI6G92ydWywAQKBQZhqgzCDSebw+vbmiSDPzturw8Y+bRvZJ0yNjcpSVEmdyOgBoGmWmAcoMUKe8qlZ//HSr/pZfKI/PUJTNoptzu+ney3oqMSbK7HgA0AhlpgHKDNDY9tJKPf3BRi3cUiZJah/n0KRR5+qGC7rIZrWYnA4A6lBmGqDMAE1btKVUT83bqG/KjkqSemfEa9rYHA1lWgQAQYAy0wBlBmherdenN74s1B8WbFP5sbrraUblpOuRMX3UtQPX0wAwD2WmAcoMcGqHj7r1xwVb9fflRfL6DDlsVt1ycTfdc2kPxUdzPQ2AtkeZaYAyA5y+rSUVemreRi3edkCSlNLOoV+N6qXrBmdyPQ2ANkWZaYAyA7SMYRhauKVUv563STsO1F1Pk9MxQdOuztFF2R1MTgcgUlBmGqDMAGfG7fHpb18W6o8Ltqqi2iNJuvK8DP2/q/oos32syekAhDvKTAOUGeDsHDrq1sy8Lfrf5UXyGVKUzaKr+nXUzUO7aUBmkiwWPn4C0PooMw1QZoDWsaW47nqaJdsP+Jed3zlR43O7aez5HRUdZTMxHYBwQ5lpgDIDtK51e45o7rJC/d+6fXJ7fJLqbrz3swsz9fMhXdUpKcbkhADCAWWmAcoMEBgHK2v01le79caXhdpXXi1JslktGpWTrvG53XRRdns+ggJwxigzDVBmgMDyeH1asKlEc5bt0pc7DvmX986I1/jcbrp2QCfFOuwmJgQQiigzDVBmgLazpbhCc/N36d1Ve3Ws1itJSoi26/rBmRqX25W7CgM4bZSZBigzQNsrr6rVP1bu1t++LFThwSpJksUi/aBXmsbndtUlPVNl5SZ8AE6CMtMAZQYwj89n6POtZZqzbJc+31rmX56VEqdxF3XVTwZ3VgLTJQBoAmWmAcoMEBx2lFXqb18W6p9f71FFTd1N+GIdNv1o4Dm6ObebeqbHm5wQQDChzDRAmQGCS2WNR++u3qvXl+3SttJK//KLe3TQ+NxuGtknnXmgAFBmGqLMAMHJMAzlf3NQc5bt0oJNJfId/9coNd6py3PSdXlOuoZ27yCnnZvxAZGIMtMAZQYIfrsPVemN5UV666siHamq9S9v57Tr+71SNSonXSN6pSkxhutrgEhBmWmAMgOEjhqPV/nfHFTexhLlbSxRaUWNf53datFF2R00qm+6RvZJ507DQJijzDRAmQFCk89naN3ecn1SUKxPNpZoe4PraySp3zmJGpWTrsv7pqtXejx3GwbCDGWmAcoMEB52lFUqb2OJPtlYolVFh9XwX68u7WPrik1OugZ3a88FxEAYoMw0QJkBwk9ZRY0+3VT3UdTi7Qf8E15KdZNeXto7TaNy0jW8Z6piHFxADIQiykwDlBkgvB2t8WjxtjJ9UlCiTzeXqvzYtxcQR0dZNbxn3QXEl/VJV/s4h4lJAbQEZaYBygwQOWq9Pn2165A+Kag7a7P3yDH/OqtFGtytvYb1SNHgrsnqn5mkOCcTYALBijLTAGUGiEyGYWjjfpe/2Gzc72q03mqR+nRM0KCuyRrUNVkDuySrc3IMFxIDQYIy0wBlBoBUdy+bhVtK9dWuw1pVeLjRWZt6afFOf7kZ1DVZfTslymG3mpAWAGWmAcoMgKbsLz+mVYVHtLLwsFYWHVbB3nJ5fI3/SXTYrerfOVEDuyZrUJdkDeyarJR2TpMSA5GFMtMAZQbA6Tjm9mr93vK6clN4SCsLD+twg7sR1+vWIVYDuyZrcNf2GtQ1WT3T2snKV8GBVkeZaYAyA+BMGIahnQeOamXhYa0qOqyVhYe1taTyhHHx0XYN6FJ/5iZJvTMSlNLOwbU3wFmizDRAmQHQWsqrarV6d901NyuLDmt10RFVub0njEuKjVLPtHbqmR6vnmntdO7xP1PjnZQc4DRRZhqgzAAIFI/Xpy0lFXXlpvCw1uw+oqJDVfI18y9rYkx9yWmnnmnx/j/TEyg5wHdRZhqgzABoS9W1Xu0oO6ptpRXaVlKprSUV2l5aqV0HjzZbcuKj7f4zOD2On9E5N72dMhKiKTmIWJSZBigzAIJBda1XOw8c1bbSSm0rOV50SitUeLBK3mZaTrzTrh7p7erO5qTFq3tanDonx+qcpBhu+IewR5lpgDIDIJjVeLzadaBKW0sqvi06pZXadeDoCV8VbygpNkrnJMXUPZLr/uycHKNzkmJ1TnKMkmOjOKuDkHa6v7+p9QBgMqfdpl4Z8eqVEd9oudvj066DRxt9VLXjwFHtPVwlV7VHR6pqdaSqVgX7XE1uN9ZhU6cmyk594UmLd/KVcoQFzswAQAiqqK7V3iPHtPfwMf+fexo8L6uoOeU2omwWdUxsXHYyEqOV0s6plHaO4386mXUcpuHMDACEsfjoKPXOiFLvjKb/ga+u9Wp/efXxclN1QtnZX16tWq+hokNVKjpUddKfFeewKSXe6S85HY6XnNT6wnN8XYd2DsU77Xy0hTYXsDLz9NNP64MPPtCaNWvkcDh05MiRE8bcd999Wrp0qTZs2KA+ffpozZo1J4xZt26dJk6cqK+++kqpqam699579fDDDwcqNgCEhegom7JS4pSVEtfkeq/PUImrWnuPHNOew1X+klPqqtGByhodqHSrrLJGbo9PR91eHT1YpcKDJy89kuS0W084s9OhQelJiolSYoNHQkyUbHzUhbMUsDLjdrt13XXXKTc3V6+99lqz42699VYtX75c69atO2Gdy+XSqFGjNHLkSL388stav369br31ViUlJenOO+8MVHQACHs2q0WdkmLUKSlGF3Rr3+QYwzBUWePRgUp3XcGpqCs6ZZVuHaz8tvTUrzvq9qrG46v72KuJSTybE++0K6FRwbE3Kjz1paepZVE2JgFFAMvME088IUmaM2dOs2NeeOEFSVJZWVmTZeaNN96Q2+3WrFmz5HA41LdvX61Zs0YzZ86kzABAgFksFsVHRyk+OqrZMzwNHXN7jxecmhMKUP3z8mO1/kf9nZMrajyqqPG0qADVi3XYGpWbdk67Yh2243/a1c5pU6zTrjjn8b877I3HOO1q57ArzmmTnWIUsoL6mpn8/Hxdcsklcjgc/mWjR4/Ws88+q8OHDys5ObnJ19XU1Kim5tuL31yupq/0BwC0nhiHTZntY5XZPva0xrs9Prmq64qNq0HJafj3xg+Pf11ljUeSVOX2qspdd33Q2XLYrd8pQ7bjJaiuGMU5bYqOsinablW0w6Zou00xDpuio6yKibLJGWVTTNTxMceXRTd47rBZuZ4oQIK6zBQXFysrK6vRsvT0dP+65srMjBkz/GeGAADByeG/vsbZ4td6vD5VVHsalR1Xda2O1nh0tMZb96e7/k/Pt8u/8/eqGq/cXp+kunJ1yOPWoaOtvad1rJa6a5nqS46zQeGJibLJYbfKabfKYa8rPg77tw9nw+c2qxz274z/7pjvbMNpsynKbpHdalWUzRJ2papFZWbKlCl69tlnTzpm06ZN6t2791mFOltTp07Vgw8+6H/ucrmUmZlpYiIAQGuy26xKjnMoOc5x6sGn4Pb4VOX2qLLGoyq3V5U1nkalqG6dV8fcHlV7fDrm9qq61qtjtV5V1/pUXVv3vNrjPb7u22XHar3+KSx8xrdnksxmt1pkt1kUZbMqymaV3Vr/d4vsx5dF2SwNllsbjK8vRfXj65b/aEBn9eucaM7+tGTwpEmTNGHChJOOyc7OPps8jWRkZKikpKTRsvrnGRkZzb7O6XTK6Wx50wcARJ66sxcOJcWefTH6LsMwVOs1VO3xqvp40TlW27AMfVt+ajxeuT0+1Xh8cnt9cnsaPLw+1dR+u/zbMV7/+hPG14/z+E7I5fEZ8vgMVdeeuO5MDeySHBplJjU1VampqYHKcoLc3Fw98sgjqq2tVVRUlCQpLy9PvXr1avYjJgAAgoXFYpHDbpHDblVCdJQpGeoLlcfnU63XUK3XJ8/xP2u9Pnl89X835PE2GNPk+G+34zn++vplPdPbmbJ/UgCvmSkqKtKhQ4dUVFQkr9frv4dMjx491K5d3Q5v375dlZWVKi4u1rFjx/xjcnJy5HA4dOONN+qJJ57QbbfdpsmTJ2vDhg16/vnn9Yc//CFQsQEACCv+QqXw/bZWwKYzmDBhgubOnXvC8oULF2rEiBGSpBEjRujzzz8/YczOnTvVrVs3SY1vmpeSkqJ7771XkydPblEWpjMAACD0MGt2A5QZAABCz+n+/g7fc04AACAiUGYAAEBIo8wAAICQRpkBAAAhjTIDAABCGmUGAACENMoMAAAIaZQZAAAQ0igzAAAgpFFmAABASKPMAACAkBawWbODSf30Uy6Xy+QkAADgdNX/3j7VNJIRUWYqKiokSZmZmSYnAQAALVVRUaHExMRm10fErNk+n0/79u1TfHy8LBZLq23X5XIpMzNTu3fvjojZuCNpf9nX8BVJ+8u+hq9I2V/DMFRRUaFOnTrJam3+ypiIODNjtVrVuXPngG0/ISEhrP/H9F2RtL/sa/iKpP1lX8NXJOzvyc7I1OMCYAAAENIoMwAAIKRRZs6C0+nU9OnT5XQ6zY7SJiJpf9nX8BVJ+8u+hq9I299TiYgLgAEAQPjizAwAAAhplBkAABDSKDMAACCkUWYAAEBIo8ycwksvvaRu3bopOjpaQ4YM0YoVK046/h//+Id69+6t6Oho9evXTx9++GEbJT07M2bM0AUXXKD4+HilpaXp2muv1ZYtW076mjlz5shisTR6REdHt1HiM/f444+fkLt3794nfU2oHldJ6tat2wn7a7FYNHHixCbHh9Jx/eKLL3T11VerU6dOslgs+s9//tNovWEYmjZtmjp27KiYmBiNHDlS27ZtO+V2W/q+bwsn29fa2lpNnjxZ/fr1U1xcnDp16qTx48dr3759J93mmbwX2sqpju2ECRNOyH7FFVeccruhdmwlNfn+tVgseu6555rdZjAf20CgzJzE22+/rQcffFDTp0/XqlWr1L9/f40ePVqlpaVNjl+2bJl+9rOf6bbbbtPq1at17bXX6tprr9WGDRvaOHnLff7555o4caK+/PJL5eXlqba2VqNGjdLRo0dP+rqEhATt37/f/ygsLGyjxGenb9++jXIvWbKk2bGhfFwl6auvvmq0r3l5eZKk6667rtnXhMpxPXr0qPr376+XXnqpyfW//e1v9cILL+jll1/W8uXLFRcXp9GjR6u6urrZbbb0fd9WTravVVVVWrVqlR577DGtWrVK//73v7Vlyxb98Ic/POV2W/JeaEunOraSdMUVVzTK/uabb550m6F4bCU12sf9+/dr1qxZslgs+vGPf3zS7QbrsQ0IA8268MILjYkTJ/qfe71eo1OnTsaMGTOaHH/99dcbY8aMabRsyJAhxl133RXQnIFQWlpqSDI+//zzZsfMnj3bSExMbLtQrWT69OlG//79T3t8OB1XwzCMX/7yl0b37t0Nn8/X5PpQPa6SjHfffdf/3OfzGRkZGcZzzz3nX3bkyBHD6XQab775ZrPbaen73gzf3demrFixwpBkFBYWNjumpe8FszS1vzfffLNxzTXXtGg74XJsr7nmGuPSSy896ZhQObathTMzzXC73Vq5cqVGjhzpX2a1WjVy5Ejl5+c3+Zr8/PxG4yVp9OjRzY4PZuXl5ZKk9u3bn3RcZWWlunbtqszMTF1zzTUqKChoi3hnbdu2berUqZOys7P185//XEVFRc2ODafj6na79fe//1233nrrSSddDdXj2tDOnTtVXFzc6NglJiZqyJAhzR67M3nfB6vy8nJZLBYlJSWddFxL3gvBZtGiRUpLS1OvXr1099136+DBg82ODZdjW1JSog8++EC33XbbKceG8rFtKcpMMw4cOCCv16v09PRGy9PT01VcXNzka4qLi1s0Plj5fD7df//9uvjii3Xeeec1O65Xr16aNWuW3nvvPf3973+Xz+fT0KFDtWfPnjZM23JDhgzRnDlz9PHHH+svf/mLdu7cqeHDh6uioqLJ8eFyXCXpP//5j44cOaIJEyY0OyZUj+t31R+flhy7M3nfB6Pq6mpNnjxZP/vZz046CWFL3wvB5IorrtDrr7+uTz/9VM8++6w+//xzXXnllfJ6vU2OD5djO3fuXMXHx+tHP/rRSceF8rE9ExExazZaZuLEidqwYcMpP1/Nzc1Vbm6u//nQoUPVp08fvfLKK3rqqacCHfOMXXnllf6/n3/++RoyZIi6du2qd95557T+304oe+2113TllVeqU6dOzY4J1eOKOrW1tbr++utlGIb+8pe/nHRsKL8XbrjhBv/f+/Xrp/PPP1/du3fXokWLdNlll5mYLLBmzZqln//856e8KD+Uj+2Z4MxMM1JSUmSz2VRSUtJoeUlJiTIyMpp8TUZGRovGB6N77rlH8+bN08KFC9W5c+cWvTYqKkoDBgzQ9u3bA5QuMJKSknTuuec2mzscjqskFRYWasGCBbr99ttb9LpQPa71x6clx+5M3vfBpL7IFBYWKi8v76RnZZpyqvdCMMvOzlZKSkqz2UP92ErS4sWLtWXLlha/h6XQPrangzLTDIfDoUGDBunTTz/1L/P5fPr0008b/b/WhnJzcxuNl6S8vLxmxwcTwzB0zz336N1339Vnn32mrKysFm/D6/Vq/fr16tixYwASBk5lZaW++eabZnOH8nFtaPbs2UpLS9OYMWNa9LpQPa5ZWVnKyMhodOxcLpeWL1/e7LE7k/d9sKgvMtu2bdOCBQvUoUOHFm/jVO+FYLZnzx4dPHiw2eyhfGzrvfbaaxo0aJD69+/f4teG8rE9LWZfgRzM3nrrLcPpdBpz5swxNm7caNx5551GUlKSUVxcbBiGYYwbN86YMmWKf/zSpUsNu91u/O53vzM2bdpkTJ8+3YiKijLWr19v1i6ctrvvvttITEw0Fi1aZOzfv9//qKqq8o/57v4+8cQTxvz5841vvvnGWLlypXHDDTcY0dHRRkFBgRm7cNomTZpkLFq0yNi5c6exdOlSY+TIkUZKSopRWlpqGEZ4Hdd6Xq/X6NKlizF58uQT1oXyca2oqDBWr15trF692pBkzJw501i9erX/GzzPPPOMkZSUZLz33nvGunXrjGuuucbIysoyjh075t/GpZdearz44ov+56d635vlZPvqdruNH/7wh0bnzp2NNWvWNHoP19TU+Lfx3X091XvBTCfb34qKCuNXv/qVkZ+fb+zcudNYsGCBMXDgQKNnz55GdXW1fxvhcGzrlZeXG7GxscZf/vKXJrcRSsc2ECgzp/Diiy8aXbp0MRwOh3HhhRcaX375pX/d97//fePmm29uNP6dd94xzj33XMPhcBh9+/Y1PvjggzZOfGYkNfmYPXu2f8x39/f+++/3/7dJT083rrrqKmPVqlVtH76FfvrTnxodO3Y0HA6Hcc455xg//elPje3bt/vXh9NxrTd//nxDkrFly5YT1oXycV24cGGT/7ut3x+fz2c89thjRnp6uuF0Oo3LLrvshP8GXbt2NaZPn95o2cne92Y52b7u3Lmz2ffwwoUL/dv47r6e6r1gppPtb1VVlTFq1CgjNTXViIqKMrp27WrccccdJ5SScDi29V555RUjJibGOHLkSJPbCKVjGwgWwzCMgJ76AQAACCCumQEAACGNMgMAAEIaZQYAAIQ0ygwAAAhplBkAABDSKDMAACCkUWYAAEBIo8wAAICQRpkBAAAhjTIDAABCGmUGAACENMoMAAAIaf8fzzz3lG67YPYAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import time\n",
    "import matplotlib.pyplot as plt\n",
    "import torch\n",
    "t0=time.time()\n",
    "\n",
    "def xyz_to_rangebearing(v):\n",
    "    r=v.norm()\n",
    "    v=v/v.norm()\n",
    "    unit_x=torch.tensor([1.0,0.0,0.0])\n",
    "    cos_bear=torch.dot(unit_x,v)\n",
    "    bear=torch.acos(cos_bear)\n",
    "    return r, bear\n",
    "    \n",
    "\n",
    "def calculate_visibility(a,bearing):\n",
    "    a_vis=torch.tensor(a)\n",
    "    fov=torch.pi/2\n",
    "    v=1/(1+torch.exp(-a_vis*torch.cos(bearing*torch.pi/fov)))#*(1+torch.exp(-torch.tensor(a_vis)))\n",
    "    v_max=1/(1+torch.exp(-a_vis))\n",
    "    v_min=1/(1+torch.exp(a_vis))\n",
    "    \n",
    "    v_normalized=(v-v_min)/(v_max-v_min)\n",
    "    return v_normalized\n",
    "\n",
    "v=torch.tensor([1.,0.,0.1])\n",
    "\n",
    "xyz_to_rangebearing(v)\n",
    "\n",
    "\n",
    "s=torch.tensor([0.0,0.0],requires_grad=True)\n",
    "pts_size=128\n",
    "pts_W=torch.rand([3,pts_size])\n",
    "\n",
    "def loss_function(s):\n",
    "    R=euler_to_rotation_matrix_zyx(angle_z=s[0],angle_y=s[1],angle_x=torch.tensor(0.))\n",
    "\n",
    "    \n",
    "    \n",
    "    pts_C=R@pts_W\n",
    "\n",
    "    loss=0.\n",
    "    for i in range(pts_size):\n",
    "        v=pts_C[:,i]\n",
    "        r,bear=xyz_to_rangebearing(v)\n",
    "        \n",
    "        if i==0:\n",
    "            vis=calculate_visibility(a=10.0,bearing=bear)*10\n",
    "            print(vis)\n",
    "        else:\n",
    "            vis=calculate_visibility(a=1.0,bearing=bear)\n",
    "        loss-=vis\n",
    "        \n",
    "    return loss\n",
    "\n",
    "\n",
    "optimizer = optim.SGD([s], lr=0.001)\n",
    "\n",
    "loss_list=[]\n",
    "s_list=[]\n",
    "\n",
    "for i in range(20):\n",
    "    optimizer.zero_grad()  # 梯度清零\n",
    "    loss = loss_function(s)\n",
    "    loss.backward()  # 反向传播计算梯度\n",
    "    #print(loss)\n",
    "    loss_list.append(loss.item())\n",
    "    s_list.append(s.detach().numpy())\n",
    "    optimizer.step()  # 更新变量\n",
    "    \n",
    "\n",
    "\n",
    "t1=time.time()\n",
    "print(t1-t0)\n",
    "print(s)\n",
    "#print(pts_W)\n",
    "R=euler_to_rotation_matrix_zyx(angle_z=s[0],angle_y=s[1],angle_x=torch.tensor(0.))\n",
    "\n",
    "pts_C=R@pts_W\n",
    "\n",
    "#pts_C_normalized=pts_C/pts_C.norm()\n",
    "#print(pts_C_normalized)\n",
    "print(loss)\n",
    "\n",
    "\n",
    "plt.plot(loss_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
